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Abstract 

We present low-temperature dynamic properties of the quantum two- 
dimensional antiferromagnetic Heisenberg model with spin S = 1/2. The 
calculation of the dynamic correlation function is performed by combining 
a projection operator formalism and the modified spin- wave theory (MSW), 
which gives a gap in the dispersion relation for finite temperatures. The so 
calculated dynamic correlation function shows a double peak structure. We 
also obtain the spin-wave damping and compare our results to some experi- 
mental data and theoretical results obtained by other authors using different 
approaches. 
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I. INTRODUCTION 



Interest in quantum antiferromagnetism is old and can be traced back to the early years 
of quantum theory, with Bethe's solution for the Heisenberg antiferromagnetic chain 1 . How- 
ever, research in this field remains very active and was further triggered by the discovery of 
high-temperature superconductivity in copper-oxide compounds. Superconductivity in the 
cuprates is attained upon doping the stoichiometric parent compounds, such as La2CuO±, 
which are believed to be experimental realizations of the two-dimensional quantum Heisen- 
berg antiferromagnet (2DQHAF) with spin S — 1/2, described by the model Hamiltonian 

w = ./E4-4 (!) 

where J > 0; denotes the nearest neighbor (NN) sites on a square lattice, without 

double counting of bonds. The 2DQHAF can be mapped into the 2D quantum nonlinear a 
model 2 (2DQNLSM) 2 and many theoretical works have been dedicated to the investigation 
of the properties of this last model. However, the mapping is rigorously valid only in the 
large S continuum limit 3 , although it can be justified on general grounds for the extreme 
quantum limit 4 S — 1/2. 

The widely held belief that antiferromagnetism plays a central role in high-temperature 
superconductivity has contributed to a noticeable proliferation of theoretical, numerical, and 
experimental works devoted to the investigation of the magnetic properties of the stoichio- 
metric parent compounds, as described by (1) and by the 2DNLSM 5 ' 6 . However, despite 
this connection with high-T c superconductivity, the understanding of the properties of the 
system is important by itself. 

Some early theoretical investigations of the 2DQHAF 7 raised doubts about the nature of 
the ground state of the model, suggesting that it would be a disordered quantum spin-liquid 
state with correlations decaying exponentially with distance 8 . Later, further investigations 
ruled out this possibility and the system is known to exhibit a broken symmetry Neel ground 
state 5 , which is destroyed by thermal fluctuations 9 when T > 0. In fact, at low tempera- 
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tures, the system is in a renormalized classical (RC) regime, that is, it behaves as a clas- 
sical 2D system with coupling constants simply renormalized by quantum fluctuations, as 
showed by Chakravarty, Halperin and Nelson 4 in their renormalization-group analysis of the 
2DQNLSM. This approach was improved by Hasenfratz and Niedermayer 10 who obtained 
an expression for the correlation length £ of the 2DQNLSM in the RC regime which agrees 
very well with experimental data for the monolayer cuprate Sr 2 Cu0 2 Cl 2 u over a certain 
temperature range. 

At higher temperatures, T/J ~ 0.5, the 2DQNLSM shows a crossover into a quantum 
critical (QC) region, with correlation length linear in 1/T 12 . Although this prediction is 
difficult to be experimentally tested in cuprates such as Sr 2 Cu0 2 Cl 2 and La 2 CuO±, given 
the high value of J in these compounds (J ~ 130meV), recent experiments on copper 
formiate tetradeuterate (CFTD) 13,14 , which is also described as a 2DQHAF 15,16 and whose 
much smaller value of J ( J ~ 6meV) allows measurements up to higher temperatures, found 
no evidence for a crossover into a QC phase. It should be remarked that the 2DQNLSM is 
expected to model the 2DQHAF in the limit of low temperatures, when the correlation length 
is very large, and the above mentioned experiments seem to settle an upper temperature 
limit for the applicability of this approach. In fact, the validity of the 2DQNLSM approach 
was shown to be inadequate to describe the behavior of the 2DQHAF with spin value greater 
than or equal to 1 in the experimentally acessible temperature region 17 . For spin S=l/2, 
a series of experiments 11 ' 18 and Monte Carlo simulations 19-21 showed that the length scale 
at which the renormalized 2DQNLSM description becomes valid is surprisingly long. The 
low energy spectrum of the S=l/2 Heisenberg model, obtained by quantum Monte Carlo 
on finite size lattices 22 , disagrees rather strongly with the prediction of the non-linear sigma 
model even when the size of the system is not too small. All these results show that it 
is important to work directly with Hamiltonian (1) instead of using the 2DQNLSM if one 
wants to make comparisons with experimental data. 

Nowadays, it is believed that the available experimental data for static properties is well 
described by a combination of low-temperature static properties 4 , QMC simulations 19 ' 20 
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at an intermediate range, and, at higher temperatures, by high-T expansion 21 and pure- 
quantum self-consistent harmonic approximation (PQSCHA) 23 . So, the attention has turned 
to the dynamical properties of the 2DQHAF at finite temperatures, since the intrinsic non- 
linearities in the equations of motion for a spin system give room for stronger quantum effects 
in the dynamics, particularly at higher excitation energies. Therefore, further microscopic 
calculations of the dynamic structure factor for (1) are still very welcome in order to allow 
a better understanding about the system's behavior at finite temperatures. 

Many of the interesting phenomena and experimental measurements in strongly corre- 
lated quantum systems are related to the dynamics of the system. Despite the considerable 
progress in many-body theory, available exact results on quantum dynamics in many-body 
systems are rather scarce. Indeed, even a systematic framework for approximate calculations 
is not well established 24 . 

In this paper, we calculate the dynamical correlation function for the 2DQHAF with 
S — 1/2 using the equation of motion approach in conjunction with projection operator 
methods, following a procedure proposed originally by Reiter 25 and further developed by 
other authors 26 . This method has proven successful in the study of the classical and quantum 
Heisenberg models in one 27 and two 28 dimensions showing good agreement with experimental 
data, Molecular Dynamic simulations and, also, with other theories. 

The calculation of the memory function, which plays a central role in the formalism, 
does not require long-range order to be valid because it depends only on correlations be- 
tween nearest neighbors. The frequency of the local spin-wave modes is also used as input 
in the method. The calculation of the static correlations and spin-wave frequency needed 
in the formalism we are using, is done within the context of the modified spin-wave (MSW) 
theory proposed by Takahashi 29 and Hirsch and Tang 30 . The MSW theory applied to an- 
tiferromagnets imposes a constraint on the total staggered magnetization to be zero in an 
isotropic system, as required by Mermin- Wagner theorem 9 . It has been shown 5 that the 
results obtained via the MSW approximation are identical to the ones obtained by Arovas 
and Auerbach 32 using a path integral formulation in the Schwinger boson representation. 
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Also, the correlation length, as calculated within MSW theory, agrees with one-loop order 
in renormalization-group theory for the 2DQNL<rM 4 . 

The combination of these two techniques, the projection operator method and the MSW 
theory was already applied by some of us 31 to study the low temperature properties of the 
quantum one-dimensional Heisenberg model with spin S = 1, where a gap is expected to 
occur. 

The dynamical structure function of (1) at finite temperature was firstly calculated by 
Auerbach and Arovas 33 as the Fourier transform of the imaginary part of the spin-density 
correlation function 

S(q,u,T) = -(S z (q,u)S z (-q,u;)), (2) 

7T 

using the Schwinger boson mean-field representation. The procedure includes processes 
where the incident particle creates quasiparticle excitations as well as scattering from ther- 
mally excited particles. However, the projection operator method, in the approximation 
proposed by Reiter 25 , goes beyond a mean field theory and, thus, the effects due to magnon 
scattering are more properly incorporated in the calculation of the dynamical structure fac- 
tor. This procedure was applied to the 2DQHAF model by Becher and Reiter 34 , using a 
standard spin-wave formalism for the calculation of the static quantities, however, as it is 
well known, conventional spin-wave theory predicts a zero gap value for finite temperatures, 
i.e., it does not take into account the inexistence of long-range order (LRO) for T > 0. 
More recently, the spin dynamics of the 2DQHAF was investigated by Nagao and Igarashi 35 
by using the self-consistent theory of Blume and Hubbard 36 . Their method is expected to 
work well at high temperatures, although the authors claim that the results obtained for 
relatively low temperatures, T/J ~ 0.4, are reliable because the relaxation function is found 
to satisfy a dynamical scaling relation consistent with the nonlinear a model analysis and, 
also, Monte Carlo simulations 37 . 

The outline of this paper is as follows: Section (II) gives a brief overview of the MSW 
results for the spin- wave energy, as obtained by Takahashi 29 . The steps leading to the 
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calculation of the dynamical structure factor are also given in this section. In Section (III), 
we discuss our numerical results and, finally, in Section (IV) we present our conclusions. 

II. THE MODIFIED SPIN- WAVE THEORY AND THE PROJECTION 

OPERATOR METHOD 

It is well known that the standard spin-wave theory is not applicable to low-dimensional 
quantum magnets at finite temperatures without modifications 9 . The consequence of the 
Mermin- Wagner theorem is enforced by hand in a variational density-matrix approach pro- 
posed by Takahashi 29 , which we will shortly review here. We start our calculations, writing 
(I) in the form 



where we divide the lattice into two sublattices A and B: spins in sublattice A (B) are 



sublattice, avoiding double counting of bonds. We now apply a Dyson- Maleev transformation 
to represent the spin operators in each antiferromagnetic sublattice in terms of bosonic 
operators 



(3) 



denoted as (Tp.) and the sum runs over all r, \ G A sublattice sites and its NN on the B 




a 




(2S - ataffja^ 




S-ala^ , 




(4) 



(5) 



following the canonical commutation relations. The Hamiltonian (3) becomes, 



H = -2NJS 2 + J £ {S (4^ + 4,^ - arte, ~ + \ a U*n - bt)%\ . (6) 



We then introduce an ideal spin-wave ansatz for the density matrix of the system 
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J2'$ indicates summation over half of the first Brillouin zone, and are given by the 
Bogoliubov transformation 

a^= cosh(0$)a$ — sinh(0£)&I- , 
Pl ? = - sinh(%)a 9 - + cosh^fcl-, (8) 

where we introduced the Fourier transform of the original boson operators, 



^-r^^- (9) 

Na + Nb = N is the number of sites in the lattice. The dispersion relation, Uq-, is obtained 
after minimizing the free energy with the constraint that the magnetization on each sublat- 
tice is zero, (Sf) = and (T2) = 0, as required by Mermin- Wagner theorem. In this way, 
we get 

^=A(1-^ 2 7?) 1/2 , (10) 
7<f = \ {cosq x + cosq y ) . 

We will write the wavevector q — (q x ,q y ) in units of the inverse lattice spacing. 

The parameters rj and A can be determined by solving the following set of self-consistent 
equations, 

2 ^ , rnj 1 J A 2 2x1/2 \ /„n 

4/ = iv^ (i-^ 7 |)v 2 2 coth |^ ( Vlci) r ( } 

The same equations have been obtained by Hirsch and Tang 30 , and, as said before, by Arovas 
and Auerbach 32 in their Schwinger boson treatment. 

Takahashi 29 was able to find out the asymptotic forms of (11) in the T — > limit and, 
also, to evaluate the 77 parameter for the spin S = 1/2 case for 4x4 and 64 x 64 lattices. 
However, for the calculation of the dynamical structure factor according to the projection 



operator procedure, we need to know the spin-wave energy for the infinite square lattice 
model at finite temperatures and, thus, we solved Eqs. (11) using an iterative numerical 
model obtaining the results displayed in Table I. 

As it is well known, there is no gap in the 2DQHAF with S — 1/2 at T — 0. Indeed, 
when T — * 0, we can see that rj — > 1. However, for finite temperatures, rj becomes smaller 
than unity and so a gap opens in the system, reflecting the fact that the correlation length, 
£, becomes finite and the long wavelength, low-energy spin-waves cannot propagate. So, 
spin-wave excitations are well defined only for wavelengths significantly smaller than £. 

Now, we describe only the main steps leading to the calculation of the dynamic structure 
factor following the projection operator method. A complete description of the theory can 
be found in Refs. 25, 26. One of the advantages of this procedure is that it allows us to 
obtain the structure factor for all values of the wavevector q, while calculations based on 
the nonlinear a model are restricted to the long wavelength limit q — > (0, 0). 

In antiferromagnets, spin-waves have two flavors, one associated with the conventional 
magnetization, M!l = S~ + and another to the staggered magnetization, VJi = 52 — T~, 
with a — x,y, z. Since the last is a lower energy mode, we will concentrate on the calculation 
of the dynamical response associated to it. We can see that the staggered magnetization is 
linear in magnon creation and annihilation operators while the uniform magnetization is a 
two-magnon operator process. In fact, the calculation of the M.^ correlation function could 
be performed by using the same procedure applied in this work, but the calculation would 
require us to go to a higher order in magnon operators. 

It is important to emphasize that we only calculate rotationally invariant quantities such 
as 71$ = | (7£;£ + 7VL + TVt). The Fourier transform of the relaxation function is given by 38 




(12) 



Here, (A, B) is the Kubo inner product of two operators A and B defined as : 



,39 




(13) 
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where < . . . > denotes the usual thermal average and (3 = l/ksT. One can show that, after 
some analytical work, the dynamical correlation function R(q, uS) is given by 

S" (to) 

R(q, u) = (7^, Ki) ( W }> f^L (14) 

[ufi-(u%+u>Vf(u>)] +[u;£» ? -H] 2 

where S^(w) and S^(cj) are the real and imaginary parts of the second order memory func- 
tion, £q-(a>), respectively. In time space, this memory function is expressed by 

^_ (QL^,e-^QL^) 
EfW " LK,) (15) 

where Q is a projection operator that projects out any term proportional to TZ^ and LIZ^, 

and L is the Liouville operator, defined by the relation LA = —i[A, H] = —iA. Reiter 25 has 

shown that, to leading order in temperature, the projection operator Q in the exponential 

function in (15) can be dropped and we can also write 

QL 2 n^L 2 TZ^-{u;l)TZ^ . (16) 

In (14), we also need to define the second frequency moment (oof), which is given by 

( W J) = {L ^ L ^\ (17) 



The second time derivatives needed to evaluate the numerator of the memory function 
(15) are directly obtained from the definition of the Liouville operator and from (3). Since 
this calculation is very straightfoward and the expressions are enormous, we will not show 
them here. We then apply the Dyson-Maleev (4) and Bogoliubov (8) transformations for 
the spin operators in those expressions. Doing so, we can simply replace the time evolution 
exp(— iLt) by the harmonic time evolution 

a^t) = e-^'o^O) (18) 
4(t)=e^4(0), (19) 

with similar equations for (3q(t) and /3t(t). So, we are left with a number of Kubo products 
of four bosonic operators which can be decoupled by means of Wick's theorem. After a 
tedious but simple calculation, we get 



S ^) = ]yE ' {A + (q,p)cos(Q + t)+A_(q,p)cos(n_t)}. 
p 



A + (q,p) and A_(q,p) are given by 

X 2 Tn t!+ n tr Je^+ - 1) 



\ 2 Tn - n - (e (3uJ '^+ - e^- ) 

where we introduced the notation q± — | ± p, and 

s(q,p) = [4J 2 ( 7(f+ + 7 <f-)(7<f + + 7<f- " 2) + (wj>] (1 + W+ ) 1/2 (l + W-) 1/2 , 
t(g,p) = [4J 2 ( 7(f+ + 7«f_)(7<f + + 7<f_ + 2) + (wj>l (1 - ^ + ) 1/2 (1 - 177,- ) 1/2 - 



(20) 



(21) 
(22) 



(23) 



In the expressions above, n,j> = (exp(/3iv$) — l) 1 is the boson occupation number and fl±(q, p) 
is defined as 



(24) 



The second moment is readily evaluated from its definition (17). It is given by the ratio 
of the following expressions 



(LftflLftg) = ' (27g + 7g + ,+ 7^) (1 _?g 2 , 1/9 coth(A(i _ ^1/2} ? (25) 



^2 2)1/2 



2T 



and 



T 1 



(26) 



A (1-777?)' 

We note that (25) is also needed in the evaluation of equation (21). 

If we take the Laplace transform of the memory function (20) and then apply the Cauchy 
formula we finally get the real and imaginary parts of the memory function as given by 



A+{q,p) + A + (q,p) + A_{q,p) + A_(g,p) 1 
u + Q + (q,p) u — Q+(q,p) u + Q-(q,p) u — Q-(q,p) J' 



and 
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A_(g,p)5(cu + + A_(q,p)5(u ~ MM)]- (28) 

Comparing equation (20) with the corresponding one obtained by Becher and Reiter (see 
equation (9) in their paper 34 ), we can note that they are very similar if we assign to 77 and 
A their zero temperature values in our expressions. But there are some slight differences 
between our definitions for s(q,p) and t(q,p), eqs. (23), and their equivalent ones due to the 
fact that, in this work, we considered two sub-lattices for the antiferro magnet. 

In the next section we discuss our numerical results for the dynamic structure factor 
calculated from the above expressions. 

III. DISCUSSION OF THE NUMERICAL RESULTS 

The evaluation of the imaginary part of the memory function is far from trivial. For this 
purpose, we adopt a numerical method introduced by Gilat and co-workers 40 and further 
developed by Wysin 41 . A typical result is shown in figure 1, for (q x ,q y ) = (ir/2,ir/2), 
t = T / J = 0.20; the frequency is given in units of JS 2 (throughout this paper, ks = h = 1)- 
A prominent feature displayed in that figure is that S^(cj) vanishes in a finite interval 
approximately centred at the zero temperature spin-wave energy. This can be understood 
if we look at eq. (28). We can see that we only have contributions when the argument of 
the delta functions, given by differences or sums of the frequency and the functions Q-(q,p) 
(which accounts for processes involving an absorption followed by the reemission of a magnon 
with a lower energy) and Q + (q,p) (that describes the decay of a magnon into two others), 
vanishes. The shape of these two functions is shown in figure 2, also for (q x , q y ) = (tt/2, 7r/2) 
and r = 0.20. We note that the non-unitary value of rj at finite temperatures, a consequence 
of the absence of long-range order, imposes the opening of a small gap between f2_ (q, p) and 
Q + (q,p), as we can see in the figure. This gap is close to the zero-temperature spin- wave 
frequency. So, as we raise the frequency, the two terms involving Q_(q,p) cease to contribute 
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to the integral in a region where the term in S(u — Q+(q,p)) still does not contribute, and 
S^(o;) is zero within it. A similar behavior was obtained by some of us in a previous work 
regarding the dynamics of the 5 = 1 antiferromagnetic chain 31 . 

The real part of the memory function is readily evaluated by a generalization 42 of a 
numerical method used in one-dimensional cases. In this way, we can directly compute 
the dynamical correlation function by means of eq. (14). In figure 3, we show the results 
obtained for (q x ,q y ) = (n/2,n/2) at r = 0.05, 0.10, 0.15 and 0.20. The results for other 
values of q are very similar, except for small wave-vectors (such as (q x ,q y ) = (tt/64, 7r/64)) 
at the highest temperatures studied. In this case, we observe that the imaginary part of 
the memory function changes its shape, from the normal one as displayed in figure 1 at the 
lowest temperature investigated (r = 0.05), into an anomalous one in higher temperatures. 
This leads to a Hurrying of the peaks observed in R(q,cu). However, this region - small 
wavevector and high temperature - is outside the validity range of the method employed in 
this work. 

We can note that the cancellation of the imaginary part of the memory function in an 
interval, as discussed above, leads to the vanishing of the dynamic structure factor in the 
same region. In particular, this feature prevents us to investigate the q — > (0, 0) limit, since 
the gap remains relatively large and there are no peaks in R(q,u). So, we cannot compare 
our results to those obtained by Becher and Reiter 34 , who found that the damping of the 
magnons is zero at T = 0. 

But, what is remarkable in the curves obtained for R(q,u) is the presence of a double- 
peak structure. It is interesting to note that the low-energy peak follows the behavior 
expected for a spin-wave peak, that is, it becomes broader as the temperature increases and 
its intensity decreases. However, the high-energy peak surprisingly has its intensity increased 
when we turn up the temperature. One can be tempted to relate this peak to two-magnons 
excitations, which should become more important when the temperature raises. It is hard to 
imagine that such a sharp structure defining a double-peak structure will be destroyed if we 
take into account higher order processes, which are expected to give important contributions 
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only at high temperatures. A similar structure was also obtained by Auerbach and Arovas 33 
in their Schwinger boson treatment of (1). Although we cannot assure that such a result, 
a double peak structure, is not an artifact of the approximations done, it should be worth 
to check if such a structure could be resolved in neutron scattering experiments. The same 
suggestion was made by Auerbach and Arovas. 

In order to calculate the damping of the spin-waves, T(q,r), we fit lorentzians to the 
data points obtained for R(q,u). T(q, r) is so obtained as the half width of the low-energy 
peak. In figure 4, we show T(q, r) as a function of the wave-vetor magnitude along two high 
symmetry directions in the antiferromagnetic Brillouin zone, for r = 0.05. We then choose 
the wave-vector (q x ,q y ) = (ir/4,ir/4), indicated in the figure, since this has an intermediate 
value for the damping, in order to compare the temperature dependence of T(q, r) with the 
experimental data obained by Thurber et a/. 43 . We show the results in figure 5. We can see 
an excelent agreement between our data, indicated by filled diamonds, and the experimental 
ones, up to a temperature of T ~ 350K. Also shown in the figure are the previous results 
by Tyc and Halperin 44 in their self-consistent calculations and we also see a good agreement 
between our results. 

IV. CONCLUSIONS 

In summary, we have calculated the dynamic structure factor for the antiferromagnetic 
Heisenberg model with S — 1/2 using a projection operator technique. This approach was 
also followed by Becher and Reiter 34 , but they used conventional spin- wave theory in order 
to calculate the static correlations needed as an input. Instead, we combined the method 
with the MSW theory, which goes beyond the linear one in the sense that it takes into 
account the inexistence of long-range order in the model at any finite temperature. We 
obtained a double-peak structure, as in a previous work of Auerbach and Arovas 33 . The 
damping, calculated by fitting lorentzians to our data points, agrees well with experimental 
data 43 and with previous theoretical calculations by Tyc and Halperin 44 up to a temperature 
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T ~ 350AT. 

It is worth to remark that, in the classical limit, the gap will vanish and the double peak 
struture will disappear. So, the double peak structure found here is a quantum effect. 

Huberman et al 45 probed the low temperature magnetic excitations of the 2D S=5/2 
AF compound Rb 2 MnF i using pulsed inelastic neutron scattering and found a dominant 
sharp peak that can be identified with one-magnon excitations. However, in addition to this 
one magnon peak, he was able to observe a relatively weak continuum scattering at higher 
energies. This continuum scattering was attributed to scattering by pairs of magnons as 
expected to happen for the M. q correlation function. This will be the subject of a future 
work. 

Acknowledgements: We acknowledge Takashi Imai for sharing the experimental data shown in 
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TABLES 



T/J 


V 


A 


0.05 


0.999999993 


2.609963279 


0.10 


0.999995178 


2.609876929 


0.15 


0.999939299 


2.599595679 


0.20 


0.999709621 


2.562049716 


0.25 


0.999280433 


2.567211935 


0.30 


0.998523449 


2.558794763 



TABLE I. Results for the rj and A parameters for some temperatures T used in our calculations. 
The values were obtained by solving Eqs (11) numerically. 
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FIGURES 
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FIG. 2. Q-(q,p) and Cl + (q,p) as a function of {p x ,Py) for {q x ,Qy) = (vr/2, 7r/2) and r = 0.20. 



20 




CO 



FIG. 3. Dynamic structure factor, R(q,oj), for (q x ,q y ) = (7r/2,7r/2) and r = 0.05,0.10,0.15 and 
20. 
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FIG. 4. Results obtained in our work (as explained in the text) for the spin-wave damping, 
T(q, r),as a function of (q 2 . + q 2 ) 1 ^ 2 for r = 0.05. Filled circles are for wavectors along the direction 
(q, 0) and open circles for (q, q). Wave-vector (q x ,q y ) = (tt/4, 7r/4) is indicated. 
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FIG. 5. Spin-wave damping, T(q, r), as a function of temperature for (q x ,Q y ) = (t/4, 7r/4) 
(filled-diamonds) . Open squares are the experimental results from Thurber et al., 43 and open 
triangles are theoretical results from Tyc and Halperin 44 . 
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